clear all
clc


data = readtable('data_boundtest');

load('estimates.mat')

J = 3;    % # of goods

id_unique = unique(data.cs);
ngoods = zeros(length(id_unique),1);

for i=1:length(id_unique)
    ngoods(i) = sum(data.cs==id_unique(i));
end

id_unique = id_unique(ngoods==3);

data = data(ismember(data.cs,id_unique),:);

N = length(id_unique);

data_wide.choice = zeros(N,1);
data_wide.x = zeros(N,J);
data_wide.z = zeros(N,J);

for i=1:N
    ind = (data.cs==id_unique(i));
    data_wide.choice(i) = find(data.chosen(ind));
    data_wide.x(i,:) = data.price(ind)';
    data_wide.z(i,:) = data.discount(ind)';
end

data_wide = struct2table(data_wide);

data_wide.x = -data_wide.x;
data_wide.z = -data_wide.z;

[data_trans] = transform_data_01_exp(data_wide,J);

s_hat = NaN(J,N);

combos_all_unique = double(combos_all_unique);

parfor i=1:N
    tic
    for j=1:J
        
        argument = [data_trans{j}.x_own(i) data_trans{j}.x_other(i,:) ...
                    data_trans{j}.z_own(i) data_trans{j}.z_other(i,:)];
        
        s_hat(j,i) = fun_sigma_est(argument,theta_NPD,combos_all_unique,m_regr,J);
        
    end
    toc
    i
end

data.s_hat = s_hat(:);

mean(data.LB<=data.s_hat&data.s_hat<=data.UB)

temp = data.s_hat-data.LB;
temp2 = temp(temp<0);
median(temp2)

temp = data.UB-data.s_hat;
temp2 = temp(temp<0);
median(temp2)

quantile(data.s_hat-data.LB,[0.05 0.15 0.25 0.5 0.75 0.95])
quantile(data.UB-data.s_hat,[0.05 0.15 0.25 0.5 0.75 0.95])

save('overid_test.mat')

ngrid = 100;
grid_q = linspace(0,1,ngrid+1);
LB_q = quantile(data.LB,grid_q);
% LB_q = [LB_q 1];
d = dataset();
d.lbquant = (1:ngrid)';
for i=1:ngrid
    temp = find(data.LB<LB_q(i+1)&data.LB>=LB_q(i));
    d.s_hat(i) = mean(data.s_hat(temp));
    d.lb(i) = mean(data.LB(temp));
    d.ub(i) = mean(data.UB(temp));
end


figure(1)
pl = plot(d.lbquant,d.s_hat','-k','linewidth',1.5);
hold on
p2 = plot(d.lbquant,d.lb','--^','color','red','markersize',4);
hold on
p3 = plot(d.lbquant,d.ub','-*','color','blue','markersize',4);
hold on
legend([pl p2 p3],'Estimated Choice Probability','Lower Bound','Upper Bound')
hold off
xlabel('Quantiles of the Lower Bound')
ylabel('Choice Probability')
ax = gca;
outerpos = ax.OuterPosition;
ti = ax.TightInset; 
left = outerpos(1) + ti(1);
bottom = outerpos(2) + ti(2);
ax_width = outerpos(3) - ti(1) - ti(3);
ax_height = outerpos(4) - ti(2) - ti(4);
ax.Position = [left bottom ax_width ax_height];
fig = gcf;
fig.PaperPositionMode = 'auto'
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
axis([min(d.lbquant)   max(d.lbquant)   0   1])




